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Modeling interacting galaxies to reproduce observed systems is still a challenge due to the extended parameter space 
(among other problems). Orbit and basic galaxy parameters can be tackled by fast simulation techniques like the restricted 
A'-body method, applied in the fundamental work by Toomre & Toomre (1972). This approach allows today for the study 
of millions of models in a short time. One difficulty for the classical restricted A'-body method is the missing orbital decay, 
not allowing for galaxy mergers. Here we present an extension of the restricted A'-body method including dynamical 
friction. This treatment has been developed by a quantitative comparison with a set of self-consistent merger simulations. 
By varying the dynamical friction (formalism, strength and direction), we selected the best-fitting parameters for a set 
of more than 250 000 simulations. We show that our treatment reliably reproduces the orbital decay and tidal features of 
merging disk galaxies for mass ratios up to ^ = 1/3 between host and satellite. We implemented this technique into our 
genetic algorithm based modeling code MiNGA and present first results. 
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1 Introduction 

Finding the initial parameters of an interacting galactic sys- 
tem is still like looking for a needle in a haystack. One chal- 
lenge is the large number of parameters describing the or- 
bital and the galactic properties. Some of these parameters 
might be derived from detailed observations, e.g. from HI 
data cubes. However, to perform an effective search in a 
high dimensional parameter space, it is necessary to use fast 
simulations and sophisticated finding strategies in parame- 
ter space. We use the code MiNGA (Theis 1999), where an 
improved restricted A'-body code is coupled to a genetic al- 
gorithm (GA). Such a strategy has also been proposed by 
Wahde (1998). For the galaxy NGC4449 Theis & Kohle 
(2001) showed, that the HI structure of a weakly interact- 
ing system can be reproduced. 

1.1 Genetic Algorithm 

MiNGA uses a genetic algorithm based on pikaia (Charbon- 
neau 1995). This kind of algorithms try to imitate nature re- 
garding the evolution of species. Heredity and mutation of 
characteristics are used to adopt the simulations to the ob- 
servations. Each model parameter is coded (normalised) to 
a gene (here we use 4 digits for a gene). All genes together 
are then forming a single string, the chromosome, which is 
fully describing a complete interaction model. The realisa- 
tion of heredity and the determination of the fitness' can be 
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done using quite different techniques. Here we use a cross- 
over operator that cuts two chromosomes at a random po- 
sition and swaps the remaining ends. Better fitting models 
are more likely parents of the next generation of models. 
This process is used to evolve the models from generation 
to generation. The fitness of the models is usually raising, 
especially if elitism^ is used. However, this evolution pro- 
cess could suffer from inbreeding^ and therefore mutation'' 
is applied. MiNGA uses either constant or changing muta- 
tion rates (depending on symbols for inbreeding). A more 
detailed description of the GA is provided in Theis (1999) 
and Theis & Kohle (2001). 

1.2 Restricted A'-body 

The first who applied the restricted A'-body method to in- 
teracting galaxies were Pfleiderer & Siedentopf (1961) and 
Toomre & Toomre (1972) - hereafter TT72. This approach 
treats the galactic centres self-consistently, while the disk 
consists of (mass-free) test particles. The main advantage of 
the restricted A'-body method is the reduction of the O(N^) 
problem of the original Newtonian equation of motion to 
about O(NNc), if Nc denotes the number of galaxies. For 
point mass galaxies the set of equations is reduced to 

r,^^^-Gj]^^.(r,-RM). (1) 

^ The best model of a generation is forwarded, if no superior model was 
found in the next generation. 

' The optimisation process got stuck caused by a too homogeneous set 
of individual models. 

^ With a low probability each chi'omosome entry might be changed. 
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r, is the position of the ;-th particle and m, its mass. Rk(t) 
describes the position of galaxy k at time t and is its dy- 
namical mass including dark matter. G represents the con- 
stant of gravitation. 

Different to TT72, MiNGA allows for a self-consistent 
description of (rigid) extended halos (Gerds 2001; Theis 
2004). Though this treatment substantially influences the 
galactic orbits (and also increases the CPU time), the re- 
stricted A^-body method, i.e. Eq. (1), can still be appUed. 

Another important process is dynamical friction. It de- 
scribes the deceleration (due to scattering) of a perturber 
moving in a background of particles. Self-consistent mod- 
eling already accounts for dynamical friction, but it is miss- 
ing in classical restricted A'-body codes and, therefore, these 
codes were not able to remodel tightly interacting or merg- 
ing systems. 

1.3 Dynamical Friction 

A simple formula for dynamical friction was derived by 
Chandrasekhar (1942) by using following assumptions: A 
point mass perturber is moving in an homogeneous, infinite 
background of particles and the mass of a background par- 
ticle is neghgible compared to the perturbers mass. 



—r- = -F{vm, cr) — In A Vm 
dt 



(2) 



The acceleration dvuldt of a massive particle M de- 
pends on the background density p, the mass of the per- 
turber M and its velocity vm- The acceleration is pointing 
opposite to the direction of the velocity, hence causing an 
effective deceleration of the particle. For more details on the 
function F{vm, cr), refer to Binney & Tremaine (1987). The 
Coulomb logarithm In A is the relation between the max- 
imum impact parameter /j^nax ™d the impact parameter 
that leads to a 90° degree deflection: 



A = 



G(M -I- m) ~ ~b^' 



(3) 



Vo denotes the velocity of the reduced particle (perturber 
and one background particle with mass m). 

Recently, efforts have been made to improve shortcom- 
ings of the approach, Eq. (2). E.g. Hashimoto et al. (2003) 
and Spirmato et al. (2003) accounted for finite halo systems. 
Just & Penarrubia (2005) focused on the influence of a den- 
sity gradient. Furthermore Jiang et al. (2008) claim, that a 
mass dependency should be applied to the Coulomb loga- 
rithm for merging time scales of dark matter galaxies in a 
cluster. 

2 Method 

We are using a set of self-consistent reference models to 
determine the appropriate formaUsm of the dynamical fric- 
tion. We have chosen isothermal spheres to serve as host 



halo galaxies. Different satellites - point masses, isother- 
mal spheres and disk galaxies - are merged with them. The 
reference models have been evolved using the gyrfalcON 
tree-code (Dehnen 2000). For an independent simulation we 
also used a direct code on a Grape6A - board (Sugimoto et 
al. 1990; Makino et al. 2003). 65 000 particles per halo have 
been used for the self-consistent models. The halo was trun- 
cated at 150kpc resulting in a mass of 5.4 • lO'^M© and a 
velocity dispersion of cr = 62km/s. The deviation between 
radial decay in our models and the reference models is used 
as a diagnostics. Detailed results of our studies will be pub- 
lished later. We varied the Coulomb logarithm, the strength 
and the direction of the dynamical friction force. The force 
itself is applied in a symmetric way to the equations of mo- 
tion of the galaxy centres. The varied parameters (C/, p and 
In A) are shown in Eq. (4). 

^ = -F{vm, o-)Cf^ {vm cos (yS) -I- sin (j8)) In A (4) 

C/ is a simple scahng factor that allows for fitting. As galax- 
ies have density gradients, the force might point not exactly 
opposite the velocity, therefore we introduce an orthogonal 
component which is adjustable via p. Finally, In A is derived 
by different approaches, these also denote our models: 

- Model A uses a constant Coulomb logarithm. 

- Model B uses a distance-dependent Coulomb logarithm 
as described by Hashimoto et al. (2003): 



lnA = ln 



(5) 



I Aba] 

tm is the distance satellite - halo centre. 

Model C uses an interpolation between two constant 

Coulomb logarithms (not presented here). 

Model D uses a mass- and distance-dependent Coulomb 

logarithm, similar to a description by Jiang et al. (2008): 



lnA = ln 



1 + 



M 



(6) 



M is the mass of the satellite and Mhaio(''M) is the mass 
of the host halo enclosed within the actual satellite's ra- 
dius ru- 

3 Results 

3.1 Isothermal satellite 

We mainly tested the merging of isothermal satellites into 
isothermal halos. In total more than 250 000 restricted A'- 
body models were compared to 20 self-consistent reference 
models (with different mass ratios). Here we present two ex- 
amples for initially circular orbits, but different mass ratios, 
i.e. q = 1/30 and q = 1/3. The first example describes a low 
mass satellite. In that case all model approaches described in 
Sec. 2 are able to reproduce the radial decay of the satellite. 
In Fig. 1 we present the best fits for models with constant 
Coulomb logarithm (model A) and for a distance-dependent 
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Fig. 1 Radial decay for an isothermal satellite within an isother- 
mal halo for a mass ratio of ^ = 1 /30: Self-consistent reference 
model generated with gyrfalcON (green solid line); best model A 
(blue dashed line) and best model B (red dash-dotted line). 
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Fig.2 Same as Fig. 1, but for a mass ratio of q = 1/3: Self- 
consistent model (green solid line); best model A (blue dashed 
line) and best model D (red dash-dotted line). 

one (model B). The latter is superior because it is able to re- 
produce the complete merging process with a deviation^ of 
6d{23Q) = 1.8 • 10"^. This example already shows the limi- 
tations of a constant Coulomb logarithm, the innermost part 
of the merging sequence occurs to quickly. 

The second example was done with a larger satellite 
mass (one third of the halo mass). In that case a constant 
Coulomb logarithm is not able to reproduce the merging 
process (Fig. 2). It either leads to a large underestimation 
of the merging time or to a different behaviour of the radial 
decay (as shown there), resulting in a deviation of 5^/(185) - 
6.8 • 10"^. Nevertheless, we were able to improve the models 
by using a mass-dependent Coulomb logarithm (model D). 
With this approach we could remodel the radial decay for 
one revolution with a deviation of ^^(frevoiution) ^ 1-5 • 10"^ 
and Sditmage) - 2.8 • 10"^ for the complete merging. 

3.2 Disk satellite 

We also caiTied out simulations using a disk-like satellite 
merging with an isothermal halo. A self-consistent model 

^ Sd(t) is derived by integrating the quadratic difference between the 
compared radial decay curves over time t. 



was set up with mkkd95 (Kuijken & Dubinski 1995) and in- 
tegrated with gyrfalcON. The details of the model param- 
eters can be found in Table 1. The determination of the 
orbital decay was done in the same manner than for the 
isothermal satellites. As we could already use our results 
from the isothermal satellites, we only needed to carry out a 
few tenth of simulations. We also compared the location of 
the disk particles, i.e. the observables. Fig. 3 shows the good 
match of the radial decay and the comparable formation of 
the trailing tidal arm. Minor mismatches like the distribu- 
tion of particles in the trailing arm might be explained by 
the different initial setup of the disk. However, the leading 
arm of the disk galaxy could not be reproduced well. This 
shows, that we need to be careful about predictions for the 
innermost regions of the merger, derived from our models. 



Table 1 Properties of the disk-merger test. 



description 


self-consistent 


improved restricted 


Host-halo: 


isoth., 65 536 part. 


isoth., static, analytic 


Disk galaxy: 


self-consistent 


restricted 


- halo: 


12 000 particles 


static, analytic 




'■halo ~ lOkpc 


'halo = 6kpc 


- disk: 


16000 part. 


16129 part. (100 rings) 


integration: 


gyrfalcON 


restricted N-body 


CPU-time: 


36min 


5.6 sec (incl. setup) 


Dyn. friction: 


self-consistent 


according equation (6) 






P = nIA, Cf = 0.89 



3.3 Genetic Algorithm run 

The last result we want to present is a complete GA run with 
the implementation of dynamical friction to improve the 
restricted A^-body code of MiNGA. The GA was provided 
with a reference model (representing a real observation), a 
merger of two disc galaxies with a mass ratio q = 1/3. Eight 
free parameters were selected - see Table 2. We have used 
our mass-dependent Coulomb logarithm of Eq. (6), where 
the strength of the dynamical friction Cf was one of the 
free parameters. The results are shown in Fig. 4. Most of 
the parameters were recovered with errors of less than 10% 
cf. also Table 2. 

4 Conclusions 

We have improved the restricted A^-body code by introduc- 
ing dynamical friction. We varied the determination of the 
Coulomb logarithm as well as the strength and direction of 
the friction force. We compared our models to self-consis- 
tent simulations in order to find the best parameterisation. 
We have shown, that radial decays of mergers up to a mass 
ratio of <7 = 1 /30 can be reliable reproduced by using a con- 
stant or distance-dependent Coulomb logarithm. With the 
introduction of more sophisticated descriptions like a mass- 
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Fig. 3 Model of a disk galaxy merging into an isothermal 
halo. Comparison between a self-consistent reference model (left) 
and an improved restricted A'-body model (right). For the self- 
consistent model only disk particles are shown (16000). The re- 
stricted model used 16 129 test particles. The initial distribution 
was set to meet optically the initial reference model. Merging was 
completed at time = 60.0 TU or 1.7Gyrs. CPU time was 36min 
for the self-consistent and 5.6 sec for the restricted model. 



Table 2 Interaction parameters that should be recovered by the 
GA. Parameter name, value for the reference model and provided 
limits are listed in columns 1 to 4. The recovered parameters can 
be found in column 5 and their relative errors in column 6. 



name 


input 


limits 


recovered 


rel. error 


Mhalo,2 


1.80 


0.54 


2.70 


1.70 


5.7 10-' 


'"disk.l 


20.00 


10.00 


30.00 


19.04 


4.8 10"^ 


'■disk,2 


3.00 


1.00 


10.00 


3.23 


7.5 10-2 


'■halo,2 


3.00 


1.00 


10.00 


3.02 


7.7 10-3 


Az 


0.00 


-2.00 


2.00 


-0.83 




Av'v 


-0.1845 


-1.00 


0.00 


-0.1500 


1.9 10-' 


AVy 


0.0537 


0.00 


1.00 


0.0650 


2.1 10-' 


Cf 


0.50 


0.10 


1.00 


0.46 


8.6 10-' 



and distance-dependent Coulomb logarithm, we were able 
to remodel radial decays for mergers up to a mass ratio of 
q = 1/3. For these models it was also essential to use an 
orientation correction of the friction force. These improve- 
ments now account for a finite system with a density gra- 
dient. However, for equal mass mergers, we were not able 
to reproduce the orbital decay. Other neglected effects like 
mass loss might be the reason for failing remodeling. A few 
recent tests including mass loss already show promising re- 
sults, so we might be able to improve the restricted A^-body 
code, again. 
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Fig. 4 Result of a Genetic Algorithm run with MINGA: Particle 
distribution of the reference model shown in the upper left panel. 
The upper right panel shows the best model after generation 1. 
Lower left: best model found by the GA after 400 generations with 
100 individuals each. Lower right: evolution of the fitness (defines 
how the original particle distribution is met) over 400 generations. 
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